{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Agricultural Pricing\n",
    "\n",
    "## Objective and Prerequisites\n",
    "\n",
    "Try this example to learn how to use mathematical optimization to tackle a common, but critical agricultural pricing problem:  Determining the prices and demand for a country’s dairy products in order to maximize total revenue derived from the sales of those products. You will learn how to model this problem as a quadratic optimization problem using the Gurobi Python API and solve it using the Gurobi Optimizer.\n",
    "\n",
    "This model is example 21 from the fifth edition of Model Building in Mathematical Programming, by H. Paul Williams on pages 276-278 and 333-335.\n",
    "\n",
    "This modeling example is at the intermediate level, where we assume that you know Python and are familiar with the Gurobi Python API. In addition, you should have some knowledge about building mathematical optimization models.\n",
    "\n",
    "\n",
    "**Download the Repository** <br /> \n",
    "You can download the repository containing this and other examples by clicking [here](https://github.com/Gurobi/modeling-examples/archive/master.zip). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem Description\n",
    "\n",
    "The government of a country wants to decide what prices should be charged for its dairy products: milk, butter, and cheese. All these products \n",
    "are produced \n",
    "(directly or indirectly) from the country’s raw milk production operations. This raw milk is divided into two main components: fat and dry matter. After subtracting the quantities of fat and dry matter, which are used for making products for export or consumption on the farms, there is a total yearly availability of 600,000 tons of fat and 750,000 tons of dry matter. This is all available for producing milk, butter and two kinds of cheese for domestic consumption. The percentage composition of the products are given in the following table:\n",
    "\n",
    "| Composition | Fat (%) | Dry matter (%) |\n",
    "| --- | --- | --- |\n",
    "| Milk | 4 | 9 |\n",
    "| Butter | 80 | 2 |\n",
    "| Cheese 1 | 35 | 30 |\n",
    "| Cheese 2 | 25 | 40 |\n",
    "\n",
    "The table below shows last year's domestic consumption (demand) and prices for the dairy products:\n",
    "\n",
    "| Dairy <br /> products | Milk | Butter | Cheese 1 | Cheese 2 |\n",
    "| --- | --- | --- | --- | --- |\n",
    "| Demand (1000 tons) | 4.82 | 0.32 |  0.21 | 0.07 |\n",
    "| Price (dollars/ton) | 297 | 720 | 1050 | 815 |\n",
    "\n",
    "The elasticities and cross-elasticities are given in the table below:\n",
    "\n",
    "| Milk | Butter | Cheese 1 | Cheese 2 | Cheese 1 to  <br /> Cheese 2 |  Cheese 2 to  <br /> Cheese 1 |\n",
    "| --- | --- | --- | --- | --- |  --- |\n",
    "| 0.4 | 2.7 | 1.1 |  0.4 | 0.1 |  0.4 |\n",
    "\n",
    "The price index cannot be raised higher than last year.This constraint establishes that the new prices must be such that the total cost of last year’s consumption would not be increased. \n",
    "Last year's price index is 1.939 (measured in thousand dollars).\n",
    "\n",
    "The goal is to determine what prices and demand  maximizes the total revenue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Formulation\n",
    "\n",
    "### Sets and Indices\n",
    "\n",
    "$d \\in \\text{Dairy}=\\{\\text{milk}, \\text{butter}, \\text{cheese1}, \\text{cheese2} \\}$\n",
    "\n",
    "$c \\in \\text{Components}=\\{\\text{fat}, \\text{dry_matter} \\}$\n",
    "\n",
    "### Parameters\n",
    "\n",
    "$\\text{capacity}_{c} \\in \\mathbb{R}^+$: Yearly availability of component $c$ (1000 tons).\n",
    "\n",
    "$\\text{qtyper}_{c,d} \\in [0,1]$: Percentage of component $c$ in dairy product $d$.\n",
    "\n",
    "$\\text{consumption}_{d} \\in \\mathbb{R}^+$: Last year domestic consumption of dairy product  $d$ (1000 tons).\n",
    "\n",
    "$\\text{price}_{d} \\in \\mathbb{R}^+$: Last year price of dairy product $d$ (dollars/1000 tons).\n",
    "\n",
    "$\\text{elasticity}_{d} \\in \\mathbb{R}^+$: Last year price elasticity of domestic consumption of dairy product $d$.\n",
    "\n",
    "$\\text{elasticity12} \\in \\mathbb{R}^+$: Last year price cross-elasticity of domestic consumption of cheese 1 and cheese 2.\n",
    "\n",
    "$\\text{elasticity21} \\in \\mathbb{R}^+$: Last year price cross-elasticity of domestic consumption of cheese 2 and cheese 1.\n",
    "\n",
    "$\\text{prcIndex} \\in \\mathbb{R}^+$: Price index reflecting last year total consumption cost.\n",
    "\n",
    "### Decision Variables\n",
    "\n",
    "$\\text{p}_{d} \\in \\mathbb{R}^+$: Price of dairy product $d$ (dollars/1000 tons).\n",
    "\n",
    "$\\text{q}_{d} \\in \\mathbb{R}^+$: Demand of dairy product $d$ (1000 tons).\n",
    "\n",
    "### Constraints\n",
    "\n",
    "**Capacity**: The limited availabilities of fat and dry matter are enforced by the following constraints.\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{d \\in \\text{Dairy}}{\\text{qtyper}_{c,d}*\\text{q}_{d} } \\leq \\text{capacity}_{c} \\quad \\forall c \\in \\text{Components}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "**Price index**: This constraint establishes that the new prices must be such that the total cost of last year’s consumption would not be increased.\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{d \\in \\text{Dairy}}{\\text{consumption}_{d}*\\text{p}_{d} } \\leq \\text{prcIndex}\n",
    "\\end{equation}\n",
    "\n",
    "**Elasticities**: The demand variables $q_{d}$ are related to the price variables $p_{d}$  through the price elasticities relationships. We approximate the elasticities with linear relationships.\n",
    "\n",
    "Milk elasticity.\n",
    "$$\n",
    "(\\text{q}_{milk} - \\text{consumption}_{milk})/\\text{consumption}_{milk}) = -\\text{elasticity}_{milk}*(\\text{p}_{milk} - \\text{price}_{milk})/\\text{price}_{milk})\n",
    "$$\n",
    "\n",
    "Butter elasticity.\n",
    "$$\n",
    "(\\text{q}_{butter} - \\text{consumption}_{butter})/\\text{consumption}_{butter}) = -\\text{elasticity}_{butter}*(\\text{p}_{butter} - \\text{price}_{butter})/\\text{price}_{butter})\n",
    "$$\n",
    "\n",
    "Cheese 1 elasticity.\n",
    "$$\n",
    "(\\text{q}_{cheese1} - \\text{consumption}_{cheese1})/\\text{consumption}_{cheese1}) = -\\text{elasticity}_{cheese1}*(\\text{p}_{cheese1} - \\text{price}_{cheese1})/\\text{price}_{cheese1}) \n",
    "$$\n",
    "\n",
    "$$\n",
    "+ elasticity12*(\\text{p}_{cheese2} - \\text{price}_{cheese2})/\\text{price}_{cheese2})\n",
    "$$\n",
    "\n",
    "Cheese 2 elasticity.\n",
    "$$\n",
    "(\\text{q}_{cheese2} - \\text{consumption}_{cheese2})/\\text{consumption}_{cheese2}) = -\\text{elasticity}_{cheese2}*(\\text{p}_{cheese2} - \\text{price}_{cheese2})/\\text{price}_{cheese2}) \n",
    "$$\n",
    "\n",
    "$$\n",
    "+ elasticity21*(\\text{p}_{cheese1} - \\text{price}_{cheese1})/\\text{price}_{cheese1}) \n",
    "$$\n",
    "\n",
    "### Objective Function\n",
    "\n",
    "**Revenue**: The objective is to maximize total revenue.\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{Maximize} \\quad \\sum_{d \\in \\text{Dairy}}{\\text{q}_{d}*\\text{p}_{d} }\n",
    "\\end{equation}\n",
    "\n",
    "---\n",
    "## Python Implementation\n",
    "\n",
    "We import the Gurobi Python Module and other Python libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install gurobipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "import gurobipy as gp\n",
    "from gurobipy import GRB\n",
    "\n",
    "# tested with Python 3.11 & Gurobi 11.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Input data\n",
    "\n",
    "We define all the input data for the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# List of dairy products.\n",
    "\n",
    "dairy = ['milk', 'butter', 'cheese1', 'cheese2']\n",
    "\n",
    "components = ['fat', 'dryMatter']\n",
    "\n",
    "\n",
    "# Create a dictionary to capture the percentage composition of the products.\n",
    "\n",
    "cd, qtyper = gp.multidict({\n",
    "    ('fat','milk'): 0.04,\n",
    "    ('fat','butter'): 0.8,\n",
    "    ('fat','cheese1'): 0.35,\n",
    "    ('fat','cheese2'): 0.25,\n",
    "    ('dryMatter','milk'): 0.09,\n",
    "    ('dryMatter','butter'): 0.02,\n",
    "    ('dryMatter','cheese1'): 0.3,\n",
    "    ('dryMatter','cheese2'): 0.4\n",
    "})\n",
    "\n",
    "# Create a dictionary to capture the yearly availability of components (1000 tons).\n",
    "\n",
    "components, capacity = gp.multidict({\n",
    "    ('fat'): 600,\n",
    "    ('dryMatter'): 750\n",
    "})\n",
    "\n",
    "# Create a dictionary to capture last year's domestic consumption and prices\n",
    "\n",
    "dairy, consumption, price, elasticity = gp.multidict({\n",
    "    ('milk'): [4.82, 0.297, 0.4],\n",
    "    ('butter'): [0.32, 0.72, 2.7],\n",
    "    ('cheese1'): [0.21, 1.05, 1.1],\n",
    "    ('cheese2'): [0.07, 0.815, 0.4]\n",
    "})\n",
    "\n",
    "elasticity12 = 0.1\n",
    "elasticity21 = 0.4\n",
    "\n",
    "priceIndex = 1.939"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Deployment\n",
    "\n",
    "We create a model and the variables. The decision variables of this model are the prices and demand for the dairy products. \n",
    "\n",
    "Solving bilinear problems with Gurobi is as easy as configuring the global parameter `nonConvex`, and setting this parameter to the value of 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using license file c:\\gurobi\\gurobi.lic\n",
      "Changed value of parameter nonConvex to 2\n",
      "   Prev: -1  Min: -1  Max: 2  Default: -1\n"
     ]
    }
   ],
   "source": [
    "model = gp.Model('AgriculturalPricing')\n",
    "\n",
    "# Set global parameters. \n",
    "model.params.nonConvex = 2\n",
    "\n",
    "# Quantity of dairy products.\n",
    "qvar = model.addVars(dairy, name=\"qvar\")\n",
    "\n",
    "# Price of dairy products.\n",
    "pvar = model.addVars(dairy, name=\"pvar\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The limited availabilities of fat and dry matter are enforced by the following constraints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Capacity constraint.\n",
    "\n",
    "fatCap = model.addConstrs( (gp.quicksum(qtyper[c,d]*qvar[d] for d in dairy) <= capacity[c] for c in components  ), \n",
    "                          name='fatCap')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This constraint ensures that the new prices must be such that the total cost of last year’s consumption would not be increased."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Price index constraint.\n",
    "\n",
    "priceIndex = model.addConstr( (gp.quicksum(consumption[d]*pvar[d] for d in dairy) <= priceIndex ), name='priceIndex')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The demand variables are related to the price variables  through the price elasticities relationships. We approximate the elasticities with linear relationships."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Elasticity constraints\n",
    "\n",
    "elasMilk = model.addConstr( (qvar['milk']-consumption['milk'])/consumption['milk']  \n",
    "                           == -elasticity['milk']*(pvar['milk']-price['milk'])/price['milk'], name='elasMilk')\n",
    "\n",
    "elasButter = model.addConstr( (qvar['butter']-consumption['butter'])/consumption['butter']  \n",
    "                           == -elasticity['butter']*(pvar['butter']-price['butter'])/price['butter'], name='elasButter')\n",
    "\n",
    "elasCheese1 = model.addConstr( (qvar['cheese1']-consumption['cheese1'])/consumption['cheese1']  \n",
    "                           == -elasticity['cheese1']*(pvar['cheese1']-price['cheese1'])/price['cheese1']\n",
    "                              +elasticity12*(pvar['cheese2']-price['cheese2'])/price['cheese2'] , name='elasCheese1')\n",
    "\n",
    "elasCheese2 = model.addConstr( (qvar['cheese2']-consumption['cheese2'])/consumption['cheese2']  \n",
    "                           == -elasticity['cheese2']*(pvar['cheese2']-price['cheese2'])/price['cheese2']\n",
    "                              +elasticity21*(pvar['cheese1']-price['cheese1'])/price['cheese1'] , name='elasCheese2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The objective function is to maximize revenue."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Quadratic objective function.\n",
    "\n",
    "obj = gp.quicksum(qvar[d]*pvar[d] for d in dairy)\n",
    "\n",
    "model.setObjective(obj, GRB.MAXIMIZE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gurobi Optimizer version 9.1.0 build v9.1.0rc0 (win64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 7 rows, 8 columns and 22 nonzeros\n",
      "Model fingerprint: 0x4fffdf2e\n",
      "Model has 4 quadratic objective terms\n",
      "Coefficient statistics:\n",
      "  Matrix range     [2e-02, 1e+01]\n",
      "  Objective range  [0e+00, 0e+00]\n",
      "  QObjective range [2e+00, 2e+00]\n",
      "  Bounds range     [0e+00, 0e+00]\n",
      "  RHS range        [1e+00, 8e+02]\n",
      "Presolve removed 2 rows and 0 columns\n",
      "\n",
      "Continuous model is non-convex -- solving as a MIP.\n",
      "\n",
      "Presolve removed 2 rows and 0 columns\n",
      "Presolve time: 0.00s\n",
      "Presolved: 14 rows, 13 columns, 36 nonzeros\n",
      "Presolved model has 4 bilinear constraint(s)\n",
      "Variable types: 13 continuous, 0 integer (0 binary)\n",
      "\n",
      "Root relaxation: objective 2.791205e+00, 11 iterations, 0.00 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0    2.79121    0    4          -    2.79121      -     -    0s\n",
      "H    0     0                       2.0422947    2.79121  36.7%     -    0s\n",
      "     0     0    2.23682    0    4    2.04229    2.23682  9.53%     -    0s\n",
      "     0     0    2.18215    0    4    2.04229    2.18215  6.85%     -    0s\n",
      "H    0     0                       2.0596849    2.18215  5.95%     -    0s\n",
      "     0     0    2.17036    0    4    2.05968    2.17036  5.37%     -    0s\n",
      "     0     0    2.16917    0    4    2.05968    2.16917  5.32%     -    0s\n",
      "     0     0    2.16441    0    4    2.05968    2.16441  5.08%     -    0s\n",
      "     0     0    2.16252    0    4    2.05968    2.16252  4.99%     -    0s\n",
      "     0     0    2.15497    0    4    2.05968    2.15497  4.63%     -    0s\n",
      "H    0     0                       2.0598300    2.15497  4.62%     -    0s\n",
      "     0     0    2.15391    0    4    2.05983    2.15391  4.57%     -    0s\n",
      "H    0     0                       2.0626181    2.15391  4.43%     -    0s\n",
      "H    0     0                       2.0653760    2.15391  4.29%     -    0s\n",
      "H    0     0                       2.0656667    2.15391  4.27%     -    0s\n",
      "     0     2    2.15391    0    4    2.06567    2.15391  4.27%     -    0s\n",
      "*   85    35              18       2.0658550    2.07220  0.31%   1.6    0s\n",
      "*  103    33              16       2.0658649    2.06832  0.12%   1.6    0s\n",
      "*  139    49              18       2.0663278    2.06781  0.07%   1.5    0s\n",
      "*  161    53              19       2.0663286    2.06781  0.07%   1.4    0s\n",
      "*  164    53              18       2.0663408    2.06781  0.07%   1.4    0s\n",
      "*  165    53              18       2.0663409    2.06781  0.07%   1.4    0s\n",
      "*  209    59              21       2.0663474    2.06758  0.06%   1.4    0s\n",
      "*  210    59              21       2.0663482    2.06758  0.06%   1.4    0s\n",
      "*  211    59              20       2.0663508    2.06758  0.06%   1.4    0s\n",
      "*  226    53              27       2.0663906    2.06712  0.04%   1.3    0s\n",
      "*  339    58              24       2.0663946    2.06654  0.01%   1.2    0s\n",
      "*  341    58              25       2.0663954    2.06654  0.01%   1.2    0s\n",
      "*  368    58              22       2.0663973    2.06649  0.00%   1.2    0s\n",
      "*  371    58              23       2.0663990    2.06649  0.00%   1.2    0s\n",
      "\n",
      "Cutting planes:\n",
      "  RLT: 9\n",
      "\n",
      "Explored 397 nodes (484 simplex iterations) in 0.11 seconds\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 10: 2.0664 2.0664 2.0664 ... 2.06634\n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 2.066398260370e+00, best bound 2.066493116553e+00, gap 0.0046%\n"
     ]
    }
   ],
   "source": [
    "# Verify model formulation\n",
    "\n",
    "model.write('AgriculturalPricing.lp')\n",
    "\n",
    "# Run optimization engine\n",
    "\n",
    "model.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "\n",
    "The table below shows the price (dollars/ton) and demand (tons) at equilibrium, for each dairy product. The total revenue generated is $\\$ 2,066,398,260$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Products</th>\n",
       "      <th>Price</th>\n",
       "      <th>Demand</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <td>milk</td>\n",
       "      <td>$321.00</td>\n",
       "      <td>4,661,067.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <td>butter</td>\n",
       "      <td>$422.00</td>\n",
       "      <td>677,334.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <td>cheese1</td>\n",
       "      <td>$839.00</td>\n",
       "      <td>264,129.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <td>cheese2</td>\n",
       "      <td>$1,116.00</td>\n",
       "      <td>54,042.00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       " Products      Price        Demand\n",
       "     milk    $321.00  4,661,067.00\n",
       "   butter    $422.00    677,334.00\n",
       "  cheese1    $839.00    264,129.00\n",
       "  cheese2  $1,116.00     54,042.00"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Output Report\n",
    "price_demand = pd.DataFrame(\n",
    "    {\n",
    "        \"Products\": dairy,\n",
    "        \"Price\": [f\"{round(1000*pvar[d].x):.2f}\" for d in dairy],\n",
    "        \"Demand\" : [f\"{round(1e6*qvar[d].x):.2f}\" for d in dairy],\n",
    "    },\n",
    ")\n",
    "price_demand.index=[''] * len(price_demand)\n",
    "price_demand"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## References\n",
    "\n",
    "H. Paul Williams, Model Building in Mathematical Programming, fifth edition.\n",
    "\n",
    "Copyright © 2020 Gurobi Optimization, LLC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
